On the horseshoe drag of a low- mass planet. I - Migration in 

isothermal disks 



J. Casoli 1 

Laboratoire AIM, CEA/DSM - CNRS - Universite Paris Diderot, Irfu/Service 
d'Astrophysique, Bat. 709, CEA/Saclay, 91191 Gif-sur-Yvette, France 

jules . casoli@cea.fr 

and 
F. S. Masset 2 

Laboratoire AIM, CEA/DSM - CNRS - Universite Paris Diderot, Irfu/Service 
d'Astrophysique, Bat. 709, CEA/Saclay, 91191 Gif-sur-Yvette, France 

frederic.masset@cea.fr 

Received ; accepted 



x Send offprint requests to jules.casoli@cea.fr 

2 Also at ICF-UNAM, Av. Universidad s/n, Cuernavaca, Morelos, CP. 62210, Mexico 



- 2 - 



ABSTRACT 

We investigate the unsaturated horseshoe drag exerted on a low-mass planet 
by an isothermal gaseous disk. In the globally isothermal case, we use a formal- 
ism, based on the use of a Bernoulli invariant, that takes into account pressure 
effects, and that extends the torque estimate to a region wider than the horse- 
shoe region. We find a result that is strictly identical to the standard horseshoe 
drag. This shows that the horseshoe drag accounts for the torque of the whole 
corotation region, and not only of the horseshoe region, thereby deserving to be 
called corotation torque. 

We find that evanescent waves launched downstream of the horseshoe U-turns 
by the perturbations of vortensity exert a feed-back on the upstream region, 
that render the horseshoe region asymmetric. This asymmetry scales with the 
vortensity gradient and with the disk's aspect ratio. It does not depend on the 
planetary mass, and it does not have any impact on the horseshoe drag. Since 
the horseshoe drag has a steep dependence on the width of the horseshoe region, 
we provide an adequate definition of the width that needs to be used in horseshoe 
drag estimates. 

We then consider the case of locally isothermal disks, in which the tempera- 
ture is constant in time but depends on the distance to the star. The horseshoe 
drag appears to be different from the case of a globally isothermal disk. The 
difference, which is due to the driving of vortensity in the vicinity of the planet, 
is intimately linked to the topology of the flow. We provide a descriptive inter- 
pretation of these effects, as well as a crude estimate of the dependency of the 
excess on the temperature gradient. 

Subject headings: Planetary systems: formation — planetary systems: protoplanetary 
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disks — Accretion, accretion disks — Methods: numerical — Hydrodynamics 
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Introduction 



Planetary migration has emerged in the last decades as one of the key processes th at 



shap e forming planetary systems. Since the early work of 



Goldreich fc Tremaind (11979 



1980I ). who gave the expressions of the tidal torque between a disk and a perturber at an 



isolated resonance, the tidal torque betwe en a protoplanet an d the gaseous disk has been 



investigated in great detail. In particular, 



Warcll ( 1986 



19881 ) worked out the sum of the 



torques at all Lindblad resonances (the so-called differential Lindblad torque) and found 
that forming planets should undergo a fast decay towards their central object. Apart from 
the differential Lindblad torque, the total tidal torque between a planet and a disk features 
another component, the corotation torque, which arises from material located in the vicinity 
of the orbit. Most of the early works on planetary migration focused on the Lindblad 
torque, as it was considered that the latter dominates the corotation torque, and therefore 
dictates the direction of the migration and the order of magnitude of the drift rate. 

Recent works have nevertheless shown the importance of the corotation region on the 
total torque exerted by the disk on the planet. Traditionally, it was estimated, for low mass 
planets, by the use of a linear theory. This amounts to evaluating the amplitude of the 
resonant waves at different azimuthal wave-numbers excited by the planet in the corotation 
region, and their feed back on the planet. Such an approach was adopted by numerous 



authors, such as 



Tanaka et al. 



Ward! (119891) 



(120021 ) . Recently, however, 



Korycanskv Pollack! (|19931). and, more rec ently, by 



Paardekooper fc Papaloizoul (j2009bl ) have stressed 



the fact that the corotation torque acting on a low mass planet is not correctly accounted 
for by a linear estimate, except at early stages after a planet is "turned on" in a disk. After 
a relatively short time scale, the linear estimates break down, and these authors show that 
the total torque seems to be better described by the sum of the differential Lindblad torque 
and of the so-called horseshoe drag. This component of the torque arises from particles that 
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are close to the corotation and undergo a U-turn in front of or behind the planet. By doing 
this, they exchange angular momentum with the planet, exerting a torque on it. As the 
region in which these particles librate is called the horseshoe region, the resulting torque was 
called the horseshoe drag. It is evaluated by budgeting the jum ps of angular momentum 
of the parti cles undergoing a U - turn. This was first studied by 



recently by iMasset 



(12001 



20021 ) 



Masset fc Papaloizou 



fl2003h : 



Ward (1991 ), and more 



Paardekooper fc Papaloizou 



( I2009bl ). This raises the question of the distribution of the torque density in the coorbital 
region. While a linear analysis shows that the corotation torque arises from an annulus 
of radial width ~ H around corotation, the horseshoe drag approach considers only the 
horseshoe region, which can be arbitrarily thin, and whose width formally tends to zero 
as the planetary mass does. One can therefore wonder if there is a torque exerted by 
the material close to the orbit, but outside of the horseshoe region, and more generally 
how to connect the linear analysis, which considers the launch of evanescent waves in the 
coorbital region, and the horseshoe drag, which does not contemplate them. Also, since 
the coorbital region and horseshoe region become of increasing importance in studies of 
planetary migration, it is useful to derive a number of properties of this region. 

With this in mind, we present in this paper I an approach, valid in a globally isothermal 
disk, in which we use a Bernoulli invariant to handle the streamlines of the flow both in 
the horseshoe region and outside of it, and in which we take into account the evanescent 
pressure waves that arise on the downstream sides of the horseshoe region, and which can 
act back on the upstream regions. We show that taking these waves into account does not 
modify the horseshoe drag estimate, but that they induce a rear-front asymmetry of the 
horseshoe region, that scales with the vortensity gradient. We also work out the surface 
density response in the corotation region, and we indicate a proper way of estimating the 
horseshoe region width. This is critical for horseshoe drag estimates based on streamline 
analysis, since this drag scales with the fourth power of the horseshoe width. Note that 
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in paper II, we present a derivation of the horseshoe drag in an adiabatic disk. We use a 
different invariant to derive it, but parts of the derivation are very similar to the derivation 
undertaken in the present paper. 

We also propose a basic model of the mechanisms taking place in the locally isothermal 
case, when the disk's temperature depends on the radius. Namely, there is creation of 
vortensity which takes place near the planet, and which results in a corotation torque 
excess that differs from the globally isothermal case, and which depends on the temperature 
gradient. Much like in the adiabatic case considered in paper II, this torque excess comes 
from an edge term of the horseshoe drag. Nevertheless, owing to the lack of an invariant 
along the streamlines in the locally isothermal case, we only provide a rough estimate of the 
excess, and discuss the differences between the predictions of this estimate and the results 
of numerical simulations. 

2. Basics 

2.1. Problem definition and notation 

The system we study is composed of a central star, of mass M#, and, orbiting around it, 
a gaseous disk, with an embedded low- mass planet, of mass M p , which has a fixed circular 
orbit of radius a and an angular frequency fl p . We assume the disk to be thin, and we use 
a 2D cylindrical geometry, averaging or integrating all the quantities over z. 

We place ourselves in the frame corotating with the planet, in which we assume 
the flow to have reached a steady state. Nevertheless, we explicitly discard in this work 
possible saturation effects. This amounts to assuming that the material which arrives at the 
upstream side of the horseshoe U-turns does so with the vortensity of the unperturbed disk. 
In practice, this means that we consider the flow after a time longer than the time required 
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to complete the horseshoe U-turns, but shorter than half the horseshoe libration time. In a 
real situation, the horseshoe flow achieves many librations over the disk's lifetime, and the 
horseshoe drag saturate s to a level that depends on the dissipative processes at work in the 



disk (e.g. 



Masset 



20011 ). Hereafter we restrict ourselves to the unsaturated horseshoe drag. 
The saturation properties, both in the isothermal and adiabatic cases, will be considered in 
a forthcoming work. 

The angular frequency of a fluid element in the disk at a distance r of the star is 
Q(r). The corotation radius, r c , is the radius where the gas is rotating at the same speed 
as the planet: fi(r c ) = Q p . The disk is made of an ideal gas, and is globally isothermal: its 
temperature is constant in time and space. Initially, the surface density is a power law of 
the radius: 

£ = £ c (r/r c H, (1) 

where S c is the surface density at corotation. 

We denote the height of the disk by H, and we use the aspect ratio h = H/ a. As the 
disk is globally isothermal, the sound speed c s is constant across it, and the aspect ratio is 
only a function of r: h — c s /[afl(r)}. 

We make use of the vortensity w = lu/H, where u is the vertical component of the 
flow's vorticity: w = V x v\ z . With our notation, we have, initially: w = w c (r/r c )~ V , with 
V = 3/2 - a. 

The particles feel the gravitational potential $ = + $ p + $j, where = —GM^/r 
is the star's potential, $ p the potential of the planet, softened to avoid computational 
problems: 

% = —GM p / (r 2 - 2ra cos + a 2 + e 2 ) 1/2 , (2) 
(where e is the softening length, typically 1/3 of H), and where $j is the potential's indirect 
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term: 

$i = — — cos (j) = qaQlr cos <fi, (3) 
a 2 

where q = M p /M*. 

The enthalpy is denoted rj. For a globally isothermal disk of temperature T c , we have 
77 = T c ln(£/£ c ). 

We make use of the two Oort's constants, estimated at the planet's location, A v and 
Bp. We often use their values for a Keplerian disk, which reads : A p = — 3Q P /A and 
Bp = Q p /A. The p subscript is mandatory to distinguish B p from the Bernoulli invariant, a 
quantity we define in section 13.21 

We label the different quadrants of the horseshoe region as represented in figure HJ We 
denote with a subscript F (R) what is in front (at the rear) of the planet, rotation-wise. In 
a similar manner, we use a subscript + (— ) for the part of the disk that has r > r c (r < r c ). 
Lastly, the upstream (downstream) part of a horseshoe leg is denoted with a subscript u 
(d). 



2.2. Basic equations 

The governing equations are the one of a non magnetized, non self-gravitating, ideal 
isothermal gas, in the corotating frame, so that the radial velocity v r is unchanged with 
respect to the inertial frame, while the angular velocity is = rQ — rQ p . With this 
notation the continuity equation reads: 

dtE + ~dr{Erv r ) + ~<9 (£O = 0, (4) 

while the Euler equations read, respectively in the radial and azimuthal directions: 

d t v r + v r d r v r + —d (t) v r - rVti- 2Vt v vs = — ~ <9 r $, (5) 

r p r L 
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0.98 0.99 1.00 1.01 1.02 

Radius r/a 



Fig. 1. — An overview of the horseshoe region. This is extracted from the standard simulation 
(see Hj), at T = 25 orbits. The planet is at = and r = 1. The corotation radius is 
represented by the dotted line. The grey scale represent the total vortensity, and the two 
white lines are the separatrices. The asymmetry appears clearly, the front part being wider 
than the rear one. We also show here the integration domains 1 to 4 used in the text. 
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dtj+Vrdrj + ^j = -^-d^, (6) 

where j = r 2 Q is the specific angular momentum. 



3. Horseshoe drag 

3.1. A first torque expression 

To evaluate the torque exerted on the planet by the disk, we simply budget the angular 
momentum on a well-chosen domain D, located between two azimuths 4>d and —<Pd, chosen 
far enough from the planet to assume that we can neglect $ p there. The radial extent of 
D is bound by two streamlines, which respectively start at (r c + xd\ (pu) and (r c — Xd\ 4>d)- 
We assume that xd is much larger than x s , the typical half- width of the horseshoe region. 
The torque exerted by the material inside this region onto the planet is then: 



Tdisk^pianct = J J Ed+Q rdrd(j), (7) 
which can be recast, using (jSj), as: 

Tdisk^pianct = J j £ ^- [d t j + v r d r j + y d^jj - rdrdcf), (8) 

which we note V from now on. In steady state, we have, using (HJ): 

- (-d^VjEj) - d r {rT,v r j) - d^rP) rdrdcj) (9) 

D r 



Hjv-dn— (p rPe^-dn, (10) 

dD JdD 

where dn is a vector perpendicular to the edge dD of the domain D, oriented outwards, and 
with a length equal to the length of the elementary interval of integration on the edge. To 
proceed in the treatment of this expression, we must rely on the shape of the streamlines. 
Namely, for r — r c = ±x&, the streamlines are nearly circular and there is therefore no 
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advected flux of angular momentum into the domain D through the circular boundaries. 
We also use the fact that el ■ dn — on these boundaries, so that we can ultimately write: 



r 



r c +x D 



[rP + r(n - tt p )Zj}dr 

L-Jr c -x D J p 

where the R superscript indicates that the integral has to be performed on the rear side 



11) 



of the domain D (radial boundary at <p = —<Pd < 0), while the F subscript indicates that 
the integral has to be performed on the front side of the domain D (radial boundary at 
= 4>d > 0). In equation ( fTTl) . the first term of the integrand represents the pressure torque 
exerted on the material enclosed within the domain D, while the second term represents the 
budget of angular momentum brought to this region by advection. Since the flow is steady 
in the corotating frame, the angular momentum of this domain is constant in time and the 
torque is therefore integrally transmitted to the planet. We recognize i n the second term of 



the integrand 



Massct 



2001 



of equation (fTTl) the classical horseshoe 



Masset &: Papaloizou 



2003 



drag expression (IWard 



Massct ct al 



2006a 



1991 



1992 



Paardekooper fc Papaloizou 



2009bl ). but this equation also shows the pressure contribution, which has been overlooked 



in previous analysis. 



3.2. A Bernoulli invariant 

We derive, from equations ([5]) and ([6]) in steady state, a Bernoulli invariant, constant 
along streamlines in steady state: 

v 2 r 2 nl 

B = - + v + ^-—IL, (12) 

where v 2 = v 2 + v\. 

Far from the planet (i.e. in <pD and —<Pd), the streamlines are purely azimuthal, and 
the fluid elements do not feel the planet's potential, so we get, neglecting the indirect term: 

B = l -r 2 (fl - %) 2 + Mr) - \r 2 n 2 p + rj (13) 
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Following iMasset fc Papaloizoul (120031 ). we derive: 



d r B = ruj(VL — Vt p ) , 



(14) 



where we used the radial equilibrium: <9 r 0* = rVL 2 — d r rj. 

Before going further and inject this into equation fill I) , we evaluate how the Bernoulli 
invariant is modified at a given location in the disk with respect to its value in the 
unperturbed disk. Considering a small perturbation 5r of radius, 5r] of enthalpy, and 5Q of 
angular velocity, the total variation of the Bernoulli invariant is: 



SB = ruo (Q — Q p ) 5r — d r r] 5r + 5r] 



(15) 
(16) 



where 5'r] — 5r] — d r rj 5r is the local perturbation of enthalpy, at fixed radius. 



3.3. Final torque expression 



We now evaluate the total torque, using equations ffTTj) and (TBI) , which allows to 
change the variable of integration to B. We use the subscripts (1) to (4) to refer to the 
quadrants of figure [fl and we get: 



r c +x D 



[rP + r(tt - tt p )Zj}dr 



\- J r c -x D 



:i7) 



J F 



-jdB+ / -jdB 

Boo (1) U JB C (2) ^ 



-jdB- I -jdB 

'Boo (3) w J B c (4) ^ 

rr c +xn rr c +x D 

+ / rP dr — rP dr, 

Jr c -x D (1)(2) Jrc-x D (3) (4) 

where we supposed the radial boundaries to be far enough in r c ± xd so we can assume the 
streamlines to be purely circular there, corresponding to a Bernoulli invariant B^- 
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We note with a subscript the values of the unperturbed flow. The perturbation of 
angular momentum and pressure are denoted Sj and SP, and we use the conservation of 
vortensity along a horseshoe U-turn in order to write S/cjL^ = H/u\ 2 , 3 y This leads to: 



Boo (1) U 
Be 



Boo (3) U 
r c +x D 

r c -XD (1) 



Uo + 8j?)dB - r - 

JB X (2) ^ 

(i - + 6j?)dB + - 

4 JBov (4) ^ 

r c -fas.D 



(j'o + + tj?)dB 
r(P + 5P fl ) dr + / r(P + + <5P fi ) dr 

Jr c —XD (2) 
rr c +x D rr c +x D 

- / r(P ~ + 5P F ) dr- r(P + + 5P 

Jr c -x D (3) Jr c -x D (4) 



(19) 



The integrals over rP simplify, between (1) and (3) on the one hand, and (2) and (4) 
on the other hand. We denote Aj the total jump of angular momentum of a fluid element 
between the inside and outside region: Aj (B) = i+(P) — j°_(B), assuming that before and 
after the jump, the fluid element has the specific angular momentum of the unperturbed 
flow corresponding to its value of the Bernoulli invariant. By construction Ajo(B) is always 
a positive quantity. 



Using equation ffT6]) . we obtain: 



R 



l^R 



S'ri 



(20) 



Since we have 8'r] = SP/H, and dr = &BrdB = dB/(ru(Q — Q p )), we also can recast the 
integrals over perturbed pressure using B as an integration variable. We are hence left with: 



+ 



+ 



/•Be 


E 


/Boo 


UJ 


rB c 


e 


' Bqo 




rB c 




'Boo 


(1) 



Aj {B)dB 



A 3o (B)dB 



Be ^ 

Boo(l) U 

Be £ 

Boo (3) u 



l m R 



S'r) 



4 (^ ^p. 



-dP 

i 

-dP 



r^'r/^E 



dB 



ruj{Vt — Q p 



Boo (2) - 



cBe 


E 




/Boo(2) 


a; 




r 


E 




/Boo (4) 


a; 


4 - %) 


dP 







dB 
dB 
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dB 



Boo (3) ru(Q - Q p ) 



+ 



dB 



7* - 



(21) 



If we denote with B s the value of the Bernoulli invariant of the separatrices, we immediately 
see that Aj is strictly equal to for B^ < B < B s (i.e. in the region that is not part of the 
horseshoe region, so that its particles do not undergo any jump), so the integrals over A jo 
reduces to their values between B s and B c . Another simplification arises from the fact that 
the integrals which were previously over 5P - i.e. the last four terms of equation fl2Tl) - and 
over 5j - the second, third, fifth and sixth terms of equation ( l2Ti) - cancel each other, and 
we are only left with, reverting the integral boundaries: 



B s E 
B c W 



A 3o (B)dB - 



b s E 

Be ^ 



Aj (B)dB, 



(22) 



Ward! (119911 ). This derivation shows that 



which is exactly the original expression derived by 
the whole contribution of the material in the vicinity of corotation, including the material 
that does not belong to the horseshoe region and circulates beyond the separatrices, reduces 
to a standard horseshoe drag integral that only involves the vortensity profile inside of the 
horseshoe region. 



3.4. A model of pressure waves 

Although the torque is of primary interest, it is also worthwhile to have a peek at the 
anatomy of the horseshoe region, and to stress its major properties. Our first step consists 
in evaluating the pressure response to the horseshoe dynamics, azimuthally far from the 
planet. 

Before going further, we need to know the response of the disk to any small arbitrary 
perturbation of vortensity, independent of the azimuth. Indeed, the horseshoe dynamics 
creates, in the downstream sides of the horseshoe region, stripes of perturbed vortensity, 
which trigger a density (or pressure) and velocity response. We denote these quantities 
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respectively with and dv^. We start with the perturbed rotational equilibrium (equation 
ED: 

-2fi(J^ + c 8 —— + ——— =0, (23) 



where the pressure is P = c^S. Deriving this equation with respect to x leads to: 

f) 2 r)V 

-2n p d x 6 V(f> + c 2 s ^— = 0, (24) 

where, for a given quantity £ we neglected all the terms in comparatively to the ones 
in <9 X <5£, since, as we shall see, the latter scale with 1/H, which is much larger than the 1/r 
scaling of the former. In order to eliminate the unknown Sv^ in this equation, we make use 
of the perturbed vortensity: 

= -j^ ( 25 ) 



so we can write: 



<5£ - = -S c ^, (26) 



K W 

where k = (2Vl p uj c ) 1 / 2 is the epicyclic frequency 

The general solution of Eq. (|26|) is the convolution product of its right hand side by 
the Green's kernel K(x), which is the solution of: 

8Z-^dl5Z = 8(x), (27) 



and whose expression is: 



K(x) = 2^e"W/ ff , (28) 



where we have specialized to the Keplerian case, for which H = c s /Q = c s /k. This kernel 
represents the pressure response to a singular perturbation of w at x = 0, of weight 
J w(x)dx = w fE c . We note that K{x) has a unitary weight: 

f +oo 



/+oo 
K(x)dx = 1. (29) 
-oo 
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Downstream of the flow, we obtain 5w straightforwardly by making use of the conservation 
of vortensity along a streamline, and assuming radially symmetric U-turns (so that fluid 
particle undergoing a jump at r = r c + x is mapped to r = r c — x). For the part behind the 
planet: 

5w R (x) = -2—wqV for < x < x s (30) 
= elsewhere. 

Similarly, for the part in front of the planet 

5w F (x) = -2—w V for > x > -x s (31) 

r c 

= elsewhere. 

If we assume x s -C H, this gives for the perturbed surface density (at the rear of the planet, 
for the sake of definiteness) : 

SE R (x) = K*Y, C — (32) 
= K * (-— VE C J (33) 



r, 



c 



Xs 1 2xV„ , 
h c ax 



o 



25" r c 



8Z«(x) = (34) 



c 



The equation f[33|) shows that the actual profile is obtained by convolution of the 
fictitious surface density profile (— 2x/r c )VS c with the kernel K of unitary weight, which 
amounts to spreading radially this fictitious profile without changing the total mass of the 
perturbation. One can check that this fictitious profile is the profile that one would obtain 
by attributing the vortensity perturbation to the perturbed surface density, without altering 
the rotation profile (i.e. keeping the profile of vorticity of a Keplerian disk). This helps 
understanding why the horseshoe drag in a pressure su pported disk does not differ from the 



horseshoe drag estimated by considering test particles (IWard 



199l|). 
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We immediately deduce the corresponding perturbation of enthalpy (which is therefore 
flat over the horseshoe region because of our hypothesis that x s <C H): 

T 2 V 

V = r^X (35) 

Similarly, we have: 

r 2 V 



3.5. Impact on the horseshoe width 

The value of the Bernoulli invariant at the separatrices is the same as its value at the 
stagnation point. Quite in contrast with the situation of a non-barotropic disk, for which 
we shall see in paper II that the value of the invariant is discontinuous at the stagnation 
point, here the value of the Bernoulli invariant is necessarily continuous in the vicinity of 
the stagnation point, and the separatrices of the horseshoe region all share the same value 
B s of the Bernoulli invariant. As a consequence, in order to ensure that it is associated to a 
Bernoulli invariant with value B s , a given separatrix must shift radially, by a value 5x s , in 
order to compensate for the variation of the enthalpy due to the evanescent waves. At the 
rear of the planet, we have, using equations fl35l) and (fIBj) . evaluated at x — x s : 



5B = = AA p B p x s 5xf + T c ||^ (37) 



c 



x s 8ABr c H 
Or, specializing to the case of a Keplerian disk: 



&? T " V (38) 



x s 3 

A similar result holds in front of the planet: 

6x? 2 



Sr R 2 

' --Vh. (39) 



x s 3 



-Vh. (40) 
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We deduce the estimate of the relative asymmetry of the horseshoe region: 

tA? S-r = -Vh (41) 

The horseshoe region therefore becomes asymmetric when there is a vortensity gradient. 
The relative asymmetry scales with the vortensity gradient and the disk's aspect ratio, and 
it is independent of the planetary mass. 



Numerical simulations 



We performed simulations of our system using the FARGO code (IMasset 
order to check the properties predicted by the above analysis. 



2000a 



b|), in 



4.1. Numerical Setup 

All our runs are performed starting with a reference run, and then changing only one 
parameter. The spatial resolutions of these runs is 700 zones in azimuth and 474 in radius, 
for a ring spawning radially from 0.8a to 1.2a. We stress that, even if the radial width of our 
mesh seems small, we apply damping conditions at the boundaries that ensure no reflection 
of the wake. The damping factor was tuned in order to ensure that this assumption is valid. 

We use a planet with a mass ratio to the primary of 8 x 10~ 6 (approximately 2.6 M®, 
for a central star of mass 1 M ). The initial conditions are presented in section I2TT1 For the 
surface density, we use a typical value of a = — 1 (this value is not realistic as it corresponds 
to a surface density increasing outwards, but it corresponds to a large vortensity gradient 
which exacerbates the effects that we wish to check). The surface density at corotation is 
S c = 2 x 10~ 3 , and the aspect ratio h is set to 0.05 at r = r c . 

The last important parameter is the potential softening length. We take it to be 30% 
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of the disk thickness, which is a rather low value. Changing this value can have significant 
effects, which are discussed in section 16.31 

We indeed observe that the horseshoe region is asymmetric, by comparing its front and 
rear width. We measure it at = ±1 rad., far enough not to feel the planet's potential, 
and once the steady state is reached (after about 25 orbits). As it is taken in steady state, 
when the separatrix is well defined, we measure the whole full width (hereafter FW) of 
the separatrix, rather than twice the upstream half-width (distance between the upstream 
separatrix and the corotation, hereafter TUHW). The difference between the two methods 
turned out to be about 3-4%, which is small, but about the same order than the effect 
we are looking for. In fact, we can see this as an effect of pressure. Here, contrary to our 
assumptions of section [3T4l we do not strictly have x s <^ H (x a ~ 0.01a), so the shape of 
the pressure wave plays a small role, with a larger amplitude on the downstream part of the 
U-turn. To get rid of this problem, we chose to use the FW, more appropriate here. 

4.2. Dependence on the vortensity gradient 

Our first step was to conduct simulations with a varying vortensity gradient, in order 
to quantify the front-rear asymmetry of the horseshoe region. We performed 50 runs with a 
value of V varying from —4.5 to 4.5. The results are presented in figure [2j As we measure 
the relative asymmetry, the discrepancy between the FW an the TUHW is obvious here, the 
effect being cumulative (for a positive V, the front FW is greater than the TUHW, while 
the rear FW is lower than the TUHW). The fit of the numerical results with our expression 
is very good. 
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Fig. 2. — Relative asymmetry as a function of the vortensity gradient. The small and 
large crosses are extracted from simulations, the large crosses representing the asymmetry 
resulting from the measurement of the FHW, while the small ones correspond to the TUHW 
(see text). The dashed line represents our analytical expression of equation (|4ip . The fit 
is very satisfactory, except at very large gradients. Note the relatively high values of the 
asymmetry: for realistic values of V, the difference in size between the front and rear region 
can amount to 15 % of the "unperturbed horseshoe width" (which is taken to be the mean 
of the front and rear widths; see section I5TTT) . 
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4.3. Dependence on the aspect ratio 

The other factor impacting equation fj4TT) is the aspect ratio. We checked this 
dependency for some values of h, as shown in figure [3j The agreement with the theory is 
good, but not perfect. It is presumably related to modifications in the topology of the flow, 
which induce small discrepancies between the runs. This was not a problem for the runs 
with a varying V, where the streamline topology remained essentially unchanged between 
the different runs. 

4.4. Pressure profile 

In section 13.41 we considered x s <C H, so that the perturbed surface density, spread 
over a distance H, is uniform across the horseshoe region. As we have seen, this assumption 
is not always true, and we try here to relax this assumption, modeling the perturbed surface 
density as a convolution product of the perturbed vortensity by the Green function K(x), 
properly normalized. We compare this to the perturbed surface density of our simulations 
in figure HI We stress that we have used, for Sw, a truncated triangular profile, in order not 
to take into account the material close to corotation that did not have time to perform a 
U-turn. Our simulated profile agrees well with the measured profile from the simulations, 
both in amplitude and position. The differences on amplitude may be imputed to an 
uncertainty in x s , which is not perfectly defined here, as the situation is not stricto sensu in 
steady state. 

4.5. Time evolution 

The asymmetry of the horseshoe zone takes some time to establish. When the planet 
is "turned on" in the disk, the upstream region remains unperturbed for a while. A 
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Fig. 3. — Relative asymmetry as a function of the aspect ratio. The dashed line represents our 
analytical expression. The points are extracted from simulations, which have the standard 
parameters, except for the aspect ratio. Again, there is a satisfactory agreement, albeit not 
as good as the dependency on V. 
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Fig. 4. — Radial pressure profile. This plot is obtained at an azimuth of <p = 1, for a run with 
the standard parameters of section I4TT1 We compare the theoretical profile obtained by con- 
volution of the vortensity perturbation by the evanescent waves kernel, properly weighted, 
as indicated by equation fl32|) . and the outputs from the simulations. The shapes are very 
similar in the corotation region (our model, which assumes that the response is taken suf- 
ficiently far from the planet to be invariant in 0, breaks down in the region of Lindblad 
resonances because of the presence of the wake). 
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typical time for pressure wave s to develop is the t i me ne eded for particles to perform 



Baruteau &: Masset 



(120081 ). this time is approximatively: 



a horseshoe U-turn. Following 
Tu-tum = Qp^^q^ 1 / 2 1 (ApBp) . This is in good agreement with our simulations, as shown 
in figure [51 Here, a small initial asymmetry (approx. 1/3 of the final asymmetry) already 
exists, which we did not investigate, owing to its transient behavior. 



4.6. Impact of the mass 

Since we considered the low-mass planet limit, we want to know to which extent 
our expression is true. To do so, we performed calculations varying the planet mass to 
primary mass ratio (see figure [6]). We see that our expression stand up to q ~ 2 x 10 ~ 5 



(appro x. 7 M®), significantly over the usual "low-mass" planet limit (IMasset et al. 



2006al ). Presumably this is because the condition x s < H can be marginally verified while 



equations ([31]) and (TJOD still hold. 



5. Discussion 

We discuss hereafter a number of points related to our results. 



5.1. A definition of the horseshoe width 



Ward! (119911 ). except that 



The torque expression of equation ([221 is the one derived by 
the integral boundaries are expressed in terms of values of the Bernoulli invariant, rather 
than x s . For practical purposes, one would wish to have an expression in terms of the 
horseshoe half width. The latter is ill-defined, however, owing to the asymmetry of the 
horseshoe zone. In order to provide an adequate definition of the average horseshoe width, 
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Fig. 5. — Relative asymmetry as a function of time, in orbits. This is a run with standard 
parameters, except that the aspect ratio is h = 0.08. The vertical dashed line represent 
T~u-tum- Here, the asymmetry boost starts around Tjj_ turni as expected. Note that the initial 
asymmetry is not 0, presumably due to the initial pressure gradient in the disk (disks with 
vanishing pressure gradients do not exhibit initial asymmetry). 
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Fig. 6. — Relative asymmetry as a function of the planet mass to primary mass ratio. 
The dashed line represents our analytical expression, given by equation (|4Tj) ). The other 
points are extracted from simulations. Note that each point is measured at a different date, 
^measure oc 1 as the locally steady state takes more time to reach for smaller planets. 
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we need to derive B s , getting rid of 8rf. To achieve this, although they share the same value 
B s , we explicitly sum Bf and Bf, which yields, to second order in x: 



B* +B 



R 



B c + 2A p B p {x F s ) 2 + S'rf + B C + 2A p B p {x«) 2 + 5'r] 



R\2 i r/ R 



B„ 



B c + 2A V B. 



x 



F\2 



X 



R\2 



p^p 



(42) 



In order to transform equation (1221) so as to get rid of the Bernoulli invariant, we write to the 
first order in x the factors of the integrand: Aj (x) = 4B p r c x and H/u — £ c /u; c (l + Vx/r c ). 
This yields: 



X 



x\[l- V— I dB 

B c V r c 

B a 



B s 



X\ 



x\ 1 + V— dB 



S c r* B-B c , 



P 2i? p J Bc 2A p B p 



(<) 2 + 



which is the same torque expression as the one of 



Ward! (Il99ll ). provided we take: 



(*f) 2 + (x?) 2 



1/2 



(43) 



(44) 



This gives an adequate definition of the horseshoe width to be used when one wishes to 
estimate the horseshoe drag by means of a streamline analysis. 



5.2. A direct estimate of the corotation torque 



As we have seen, taking into ac count the pr essure effects does not unveil a different 



torque value from the one derived by 



Ward! (jl99ll ). The main reason that we have mentioned 



in section 13.41 is the fact that the launch of evanescent waves in the coorbital region renders 
radially more diffuse the perturbation of surface density (which would otherwise be bound to 
the horseshoe region), b ut leaves its linear den s ity un changed. This result is not surprising, 



since the fluid model of 



Goldreich &: Tremaind (119791 ) had no pressure dependence. 
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One can wonder, however, why the horseshoe region asymmetry described in section 13.51 
does not have any impact on the torque. In order to gain some insight into the reasons for 
that, we examine here how this asymmetry alters the front and rear components of the 
torque by a direct, albeit approximate, estimate of the corotation torque. 

The horseshoe U-turns, by conserving the vortensity, produce stripes of perturbed 
vortensity in front of the planet and behind it, in a disk with a vortensity gradient. These 
perturbation of vortensity yield perturbations of surface density. Disregarding the radial 
spread of these perturbations of surface density, we can write their linear mass as: 

A = / 5E(x)dx = (45) 



in front of the planet, and: 

A = (46) 
a 

behind the planet. We evaluate the torque arising from a stripe of linear mass A onto the 
planet, assuming that it starts at a distance d from the latter: 

r°° GM n ^V<7 
r stripe « a J r^Kdl = -IL_*A. (47) 

We note that the effects of the stripe in front of the planet, and of the stripe that is located 
behind it, are cumulative, since they involve perturbations of surface density of opposite 
signs, hence the corotation torque, in total, is: 



(x F ) 2 + (x R ) 2 

r « S%a d qWc — d { s> ■ (48) 



This allows to understand why the horseshoe drag remains unchanged when the horseshoe 
region becomes asymmetric: xf increases (decreases) while decreases (increases), but 
their sum (or the sum of their square) remains constant to lower order in V. Put into 
simple words, taking for example a positive torque excess (V > 0), the widening of the 
front region produces a supplementary torque excess, since the front region is producing a 
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positive torque excess. But the torque excess produced by the rear region is also positive, 
so the shrinking of this region induces a smaller torque excess, which balance, to the first 
order, the supplementary torque excess of the front region, since \5xf\ = \Sxf\. 



As a side result, we compare equation (1481) to the expression of IWardl (If 99fl ). which 
reads with our notation, in a Keplerian disk: 

3 



(49) 



Noting that x 2 s ~ a 2 q/h ( iMasset et al 



2006a : 



ir 



Paardekooper fc Papaloizoull2009al ). and 



assuming that the stripes originate at a distance ~ H from the planet, equation (I48j) can 
be recast as: 

(50) 



This is the same expression except for a 3/8 factor, presumably originating from our simple 
estimates of d and x q . 



6. Relaxation of the globally isothermal disk hypothesis 

In a globally isothermal disk, the asymmetry of the horseshoe region does 
not involve additional generation of perturbed surface density, because there is no 
generation of vortensity, since the gradients of pressure and density are always aligned. 
However, in a locally isothermal disk, where the temperature depends on the radius 
(T = T(r) = T c (r/r c )~ T ), this is not true anymore, and the vortensity created corresponds 
to an additional perturbation of surface density that exerts a torque onto the planet. For 
the sake of definiteness we call this torque the non-isothermal torque excess. For a proper 
formulation of this torque, one could follow the method of paper II, using an adequate 
invariant along the streamlines, whenever it exists, in order to derive the downstream 
vortensity and perturbed surface density. In the present situation, nevertheless, the 
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variation of the state variables of a fluid element that goes from infinity to the stagnation 
point depends on its path, so we cannot use the method of paper II for adiabatic flows. 
We expect therefore the non- isothermal torque excess to display a strong sensitivity to the 
topology of the streamlines. As we shall see hereafter, this is indeed the case. 

Prior to a fine tuned description of the excess inferred from numerical simulations, we 
give hereafter an oversimple estimate of the magnitude of the offset, which should apply 
to situations in which there is only one stagnation point, and which is adapted from a 
discussion provided in paper II, in which we interpret the adiabatic torque excess as arising 
from a constitutive asymmetry of the horseshoe region. 



6.1. An oversimple approach 

We know from our simulations that a locally isothermal disk exhibits a non-isothermal 
torque excess, even for a vanishing vortensity gradient. In order to estimate this excess, we 
provide hereafter an estimate of the intrinsic asymmetry of the horseshoe region of a locally 
isothermal disk. All this study is done for V = 0, in order to neglect the asymmetry of the 
horseshoe region acquired on the long term under the feed back of the evanescent waves, 
which does not plays a role in the torque, as explained before. 

To estimate the intrinsic asymmetry, we consider a fluid element belonging to a 
separatrix, that starts at (x s , <f>oo). By definition the torque applied to it is exactly sufficient 
to bring it to the stagnation point. For the globally isothermal case, in front of the planet, 
this reads: 



Ajgi = / r G i(t) dt 

Jto 



(51) 
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r GI (0) (dt/d<f>) d<b, (52) 

where the GI index stands for "globally isothermal". We then approximate (dt/d<ft) = 

~ a/2A p x Si Gi, and since Aj = 2aB p x s> Qi, we have: 

r<t>s 

AA p B p x 2 s = r(0) # (53) 

•/ (poo 

A similar relation holds in locally isothermal case, with a modified x s = x Si gi + $ x s and 
T = T GI + 5T. The only contributions in 5T is the one of the pressure, assuming that the 
fluid element follows a similar path in both cases: 

_ <VP 

i,] ~ E 



= ST^, (54) 

GI n 



where the LI index stands for "locally isothermal", and where 5T = Tli — Tqi- In 
equation (1511) we have made use of the fact that d^T cancels out. From equation (153")) . we 
have: 

8A p B p 5xf = [^ 5T^d<j). (55) 

We simplify further by assuming "square", instantaneous U-turns: the particles go from 
(x s , 0oo) to (x 8 ,<j> s ), then to (—x 8 ,<j> 8 ), and lastly to (— x s ,0oo). This ensures that is 
constant along it, and we are left with: 

8A p B p 5x F s ~ 5T(x s ) / — (56) 

J<j>oc ^ 

= -r— T , (57) 

a L c 

where we have replaced 5T by its first order expression — tTq x s /a. Denoting with 
AP = P s — Pqo, this reduces to: 



5 F - T AP (58) 
8A p B p a S c 

On the rear of the planet, the sign of 5T changes, and we have Sxf = —6xf. The crude 
estimate presented here gives an order of magnitude of the intrinsic asymmetry of the 
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horseshoe region, that gives rise to a torque for the reasons presented in detail in paper II. 
Since the enthalpy tends to fill the planet's potential well, we always have AP > 0, and 
the sign of 5xf only depends on r. For example, for a disk with decreasing temperature 
(r > 0), 8xf is negative, and the planet feels a negative torque excess from the front region. 
The rear region being in turn broader, the planet feels a negative torque excess from the 
rear region also, hence a globally negative torque excess. To quantify this excess, we use 
the same method as in paper II, which yields: 

5T = 2S C v+SxiAj ~ 2tAPx 2 s , (59) 

where we take v$ = 2A p x s , and Aj = 4B p ax s . 

We therefore find a non-isothermal torque excess that scales with the temperature 
gradient, and which should have a stronger dependency at lower softening lengths, since the 
potential well of the planet is then deeper, and AP is therefore larger. 

6.2. Comparison to numerical simulations 

We performed two sets of 30 numerical simulations with a varying temperature slope r 
between —3 and +3 to investigate the behavior of the non-isothermal torque excess. The 
first set uses a softening length of e = 0.3, while the second one uses e = 0.5. Due to a strong 
generation of vortices when V = 0, we adopted a non-vanishing vortensity gradient, and we 
chose V = 1. The results of these simulations are presented in figured We measured the 
quantity Ar as explained in figure The value of AP was measured by averaging AP F 
and AP^. 

There is a poor fit of the numerical data with our crude estimate. While the slope has 
the expected value (within an error range of ±25%), there are discontinuities in the case 
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Fig. 7. — Measurement of AI\ For a given run (here, r = —1.35, e = 0.3), we monitor the 
total torque r tot as a function of time. The differential Lindblad torque r Lindb i ad and the 
linear co-rotation torque T CR only take a few orbits t o establish, then the horseshoe dra g 
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develops, taking place of the linear corotation torque (jPaardekooper fc Papaloizou 
We measure Ar = <5Fhs — Tor, which, as Tor does not depends on r, should trace the 
dependency of the horseshoe drag on r. The values of r L i n dbiad + Tor and r L i n dbiad + ^hs 
are measured by averaging r tot between 2 and 3 orbits on the one hand, and between 30 and 
35 on the other hand. Slight changes of theses values did not bring differences. 
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e = 0.3 for \r\ ~ 1, which seem to disappear for the larger value of e. This is linked to the 
appearance of a "trapped region" , which we explain in the next section. 



6.3. A note on the softening length 

The streamlines in the inner part of the horseshoe region, near the planet, exhibit 
different behaviors as a function of the softening length, as shown in figure We can 
distinguish two cases: 



• a case where there are three stagnation points, as discussed in lMasset et al.l (l2006af ). 
with almost the same radius r c . The O-point is associated to the planet, being located 
almost at the same azimuth. One of the X-points lies on the separatrices. We call 

it the outer X-point. Similarly, the other X-point is the one we refer to as the inner 
X-point. 

• a case where there is only one X-point, at the intersection of the separatrices. 

One can continuously go from the first situation to the last one by varying the softening 
length. The inner X-point and the O-point merge at some point, yielding a one stagnation 
point configuration. A lower softening length favors the existence of an inner X-point, but 
it is not the only parameters that determines this property of the flow. The temperature 
gradient also plays a role (a smaller gradient favors the existence of a inner X point, while 
at large gradient we have a one stagnation point situation). We have seen this for two 
values of e, and present this for varying values of r in figure [TUJ The existence of three 
stagnation points appears to be correlated with the discontinuities at r = ±1. We comment 
on this behavior in the next section. 
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Fig. 8. — Non-isothermal torque excess. We plot this quantity as a function of r, for two 
sets of runs: e = 0.3 (X), and e = 0.5 (+). The dashed line represents our crude estimate 
of equation (J59l) . A proper linear regression fit of the simulated points would yield a slope 
within ±25% of our prediction. The break of continuity that we see for e = 0.3 exists also 
for e = 0.5, but is much less apparent, as it takes place at r w ±0.5, and results in a much 
smaller jump of the corrected torque excess. 
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6.4. The different regimes in the locally isothermal case 

The production of vortensity takes place near the planet, during the horseshoe U-turns, 
where the gradients of pressure (dominated by the planet's potential) and of temperature 
are not aligned. For the sake of definiteness, we assume that there is a negative temperature 
gradient. Therefore a negative amount of vortensity is produced behind the planet and a 
positive amount of vortensity is produced in front of it. In an idealized situation with a 
unique X-point centered on the planet, the negative perturbed vortensity is fully advected 
to the rear of the planet, while the positive perturbed vortensity perturbation is advected 
towards positive azimuth. They respectively yield a positive surface density perturbation 
behind the planet (for the reasons exposed in section 13.41 which are valid azimuthally far 
from the planet, where the gradients of temperature and density are aligned) and a negative 
perturbation of surface density in front of the planet (hence a negative non-isothermal 
torque excess). 

In a realistic situation, however, the fact that the X-point is not centered on the planet 
changes this picture. We take the example of r ~ — 1 (one can refer to figure [9b to see 
a streamline map that looks similar to this situation). The material at the rear of the X 
point is advected to the rear of the planet, as expected. But some of the material (and the 
resulting underlying created positive vortensity) at the front of the X point, but still at 
the rear of the planet, is eventually advected to the front of the planet. This lowers the 
torque excess exerted on the planet, as a part of the positive vortensity created at the rear 
of the planet and advected to the front cancels out with the negative vortensity created 
and advected to the front. This explains why the simulated points in figure M are below our 
prediction for r < — 1. 

As \t\ decreases, a "trapped region" appears, enclosed between the two X-points, and 
centered on the O-point. The material in this region orbits around the planet, and does 
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not contribute to the torque. In fact, the created vortensity supported by this material is 
trapped here, and even eventually flattens out (under the action of phase mixing, which 
requires very long times). For r < 0, the azimuth of the O-point is positive, and the 
material trapped in this region essentially has a negative created vortensity. Hence, when 
the trapped region appears, a significant part of the advected negative vortensity disappears 
from the torque, resulting in a smaller torque excess. This likely explains the discontinuity 
in figure [8] at r ~ — 1 . 

Similar arguments hold for the case r > (we then have to change the sign of the 
vortensity produced, as well as the sign of the azimuth of the outer and inner X-points, as 
can be seen in figure [10]) . 

For \t\ > 1, Ar decreases with r as announced, with a slope that is roughly the one 
predicted in section [6TTT . For |t| < 1, however, this slope is reverted. The reason for this 
behavior is unclear, and its study requires an in-depth streamline analysis, as well as a 
quantitative estimate of the vortensity produced in the vicinity of the separatrices, which is 
beyond the scope of this qualitative description. 

7. Conclusion 

We have derived an expression for the horseshoe drag that takes into account the effects 
of pressure. It yields a result identical to previous estimates of the horseshoe drag, except 
that it explicitly takes into account the material beyond the horseshoe separatrices. The 
horseshoe drag therefore accounts for the torque arising from the whole corotation region, 
and therefore should deserve to be called corotation torque. We find that the horseshoe 
region is asymmetric, in a manner that depends on the vortensity gradient. Specifically, the 
part in front of the planet is wider for a decreasing vortensity profile in the disk, while the 
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part behind the planet is more narrow. This asymmetry can be interpreted as an effect of 
evanescent pressure waves, launched from the downstream part of the horseshoe U-turns, 
which perturb the upstream region. Using a Bernoulli invariant, we derived this asymmetry 
for a globally isothermal disk, and found it to scale with the vortensity gradient and the 
disk's aspect ratio. 

The impact of the evanescent waves on the horseshoe width can easily be interpreted as 
a consequence of an azimuthal pressure gradient, whose torque counteracts or adds to the 
torque required to perform a U-turn, resulting either in a wider or in a narrower horseshoe 
region. 

This asymmetry does not alter the standard horseshoe drag. Basically, the 
supplementary torque excess resulting from the widening of one side of the horseshoe region 
is balanced by the torque deficit resulting of the shrinking of the other side of the horseshoe 
region. 

As the concept of horseshoe drag h as proved of increasi n g importance for planetary 



migration in a number of recent works (IMasset et al. 
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20081 ). it is crucial to adopt a proper estimate of the correct horseshoe width when evaluating 
the corresponding horseshoe drag. We find that the correct mean width of the horseshoe 
region should be estimated as given by equation fj44|) . As the asymmetry is relatively mild 
for typical vortensity gradients and aspect ratios, a simple arithmetic mean of the front and 
rear widths also yields estimates that are sufficiently accurate. 

Our analysis is based on the use of a Bernoulli invariant, which is made possible by 
the hypothesis that the disk is globally isothermal, hence barotropic. Such an approach is 
not possible for a locally isothermal disk. In such a situation, vortensity is created near the 
planet, which gives rise to a supplementary torque. We have presented a qualitative study 
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of this supplementary torque, whose dependence on the temperature gradient appears to 
be intimately linked to the topology of the flow in the planet vicinity, and which exists 
in a disk even without a vortensity gradient. We also provide a crude estimate of this 
dependence, which is meant to apply to situations with only one stagnation point (i.e. at 
large temperature gradient). We stress that this supplementary torque is not strong enough 
to reverse planetary migration in a disk with a surface density that decreases outwards as a 
power law of radius, even with unrealistic temperature gradients. Nevertheless, it should be 
considered whenever accurate estimates of the migration rate are required, as its amplitude 
can amount to about one third of the total horseshoe drag. A systematic study of these 
non-isothermal effects, and their generalization to the three-dimensional case, should be 
undertaken in order to provide reliable estimates of this supplementary torque, that we 
have called the non-isothermal torque excess. 

An important feature of this excess is that it is linked to an edge effect of the horseshoe 
drag. The bulk horseshoe drag, which corresponds to the classical expression, always scales 
with the vortensity gradient. Any dependence on another parameter of the disk (such as 
the temperature gradient) manifests itself as evanescent waves excited at the downstream 
separatrices, where they yield a vortensity sheet as the most tangible imprint on the flow. 
In the case of a locally isothermal disk, the lack of an invariant and the strong dependence 
on the topology of the flow render the situation very complex and hardly tractable. In the 
case of an adiabatic flow, one can exhibit, under certain circumstances, an invariant of the 
flow. This allows a rigorous analysis of the horseshoe drag and of the main properties of the 
horseshoe region. This analysis is presented in paper II. 

The numerical simulations performed in this work have been run on the 92 core cluster 
funded by the program Origine des Planetes et de la Vie of the French Institut National 
des Sciences de VUnivers. Partial support from the COAST project {Computational 
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acknowledge partial support from CONACYT project number 24936. 
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Fig. 9. — Appearance of a trapped region. This is a close overview of the horseshoe region, 
at T=13 orbits, close to the planet, which is at located at = and r — 1. The temperature 
gradient is r = —0.9. The left panel show the e = 0.3 case, while the right panel shows the 
e = 0.5 case. The thin dotted line represents the corotation. The underlying field is the 
perturbed vortensity. The thin white lines are the streamlines, and the bold ones represents 
the separatrices. We clearly see the advection of the vortensity created, and the advection of 
negative and positive vortensity to the front region (since, here, the outer stagnation point 
is behind the planet). When the inner stagnation point disappears, the material executing a 
U-turn in front of the planet passes faster by the planet, hence a smaller creation of negative 
vortensity (hardly seen with the grey scale, but there is a difference of about one half), while 
the creation of positive vortensity behind the planet remains of the same order. 
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Fig. 10. — Azimuth of the stagnation point(s) as a function of the temperature gradient, 
at T=5 orbits. The variation of the radial position is much smaller, although it exhibits a 
reproducible behavior (the front (rear) X stagnation point being always at r < a (r > a) ). 
For \r\ < 1.5, there is one O point, and two X-points, as indicated on the figure. At larger 
times, the inner X-point and the O-point tend to merge, reducing the interval upon which 
there are three stagnation points. Nevertheless, it remains quite large (\t\ < 1) in steady 
state. 
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